function [cy1shock,cy2shock,w1shock,w2shock,Vshock] = Value_Iteration_demo(cy1,cy2,w1,w2,V,Vtilde)

global wgrid w_thr wmax wmax_shock betta delta gamma pi1 pi2 e1_shock e2_shock ey1 ey2 v1aut v2aut nw n_shock

crit=1e-10;
optionsset=optimoptions('fmincon','Display','off','TolX',10^(-8),'TolFun',crit,'Algorithm','sqp');
options=optimoptions('fsolve','Display','off');

%Preallocation Vector
cy1shock=zeros(nw,1);
cy2shock=zeros(nw,1);
w1shock=zeros(nw,1);
w2shock=zeros(nw,1);
Vshock=zeros(nw,1);

function OBJ_neg = value_unconstrained(c) %c(1): cys1, c(2): cys2, c(3): omegas1, c(4): omegas2 
    
            OBJ_neg=...
            -(pi1.*(utility(c(1),gamma)+(betta./(delta.*n_shock)).*utility(e1_shock-n_shock.*c(1),gamma)) ...
            + pi2.*(utility(c(2),gamma)+(betta./(delta.*n_shock)).*utility(e2_shock-n_shock.*c(2),gamma)) + ...
            (delta.*n_shock).*(pi1*interp1(wgrid,V,c(3),'spline',-1e10) + pi2.*interp1(wgrid,V,c(4),'spline',-1e10)));         
  end

 function [ineqcon,eqcon] = constraints_fun(c,omega)
        eqcon=zeros(1,1);
        ineqcon=zeros(2,1);
        
        eqcon(1)=v2aut-utility(c(2),gamma) - betta*c(4); %also equality
        ineqcon(2)=v1aut-utility(c(1),gamma) - betta*c(3);
        ineqcon(1)=omega-(pi1*utility(e1_shock-n_shock.*c(1),gamma)+pi2*utility(e2_shock-n_shock.*c(2),gamma));       
 end

lb = [0;0;w_thr;w_thr];
ub = [ey1;ey2;max(wmax,wmax_shock);max(wmax,wmax_shock)];  
startguess=[cy1 cy2 w1 w2];

for k=1:nw
    guess=startguess(k,:)';
    [out,Vshock(k,1)]=fmincon(@(x)value_unconstrained(x),guess,[],[],[],[],lb,ub,@(x)constraints_fun(x,wgrid(k)),optionsset);
    cy1shock(k,1)=out(1);
    cy2shock(k,1)=out(2);
    w1shock(k,1)=out(3);
    w2shock(k,1)=out(4);
    
end
  Vshock=-Vshock;
end

